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FOREWORD 


This is a progress report on the research project "Numerical Solutions 
of Three-Dimensional Navier-Stokes Equations for Closed-Bluff Bodies." 
Specific efforts were directed in the area of "A Conservative Approach for 
Flow Field Calculations on Multiple Grids." 

The period of performance on this research was July 15, 1987 through 
March 31, 1988. This work was supported by the NASA Langley Research Center 
through Cooperative Agreement NCC1-68. The Cooperative Agreement was 
monitored by Dr. Robert E. Smith, Jr. , of the Analysis and Computation 
Division (Computer Applications Branch), NASA Langley Research Center, 


MS/125. 
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A CONSERVATIVE APPROACH FOR FLOW FIELD 
CALCULATIONS ON MULTIPLE GRIDS 


By 


M. Kathong^ and S. N. Tiwari^ 


SUMMARY 


In the computation of flow fields about 
complex configurations, it is very difficult to 
construct body-fitted coordinate systems. An 
alternative approach is to use several grids at 
once, each of which is generated independently. 
This procedure is called the "multiple grids" 
or "zonal grids" approach and its applications 
are investigated in this study. The method 
follows the conservative approach and provides 
conservation of fluxes at grid interfaces. The 
Euler equations are solved numerically on such 
grids for various configurations. The 
numerical scheme used is the finite-volume 
technique with a three-stage Runge-Kutta time 
integration. The code is vectorized and 
programmed to run on the CDC VPS-32 computer. 
Some steady state solutions of the Euler 
equations are presented and discussed. 
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1. INTRODUCTION 

Basically, there are three approaches or methods which can be used to 
solve a problem in fluid mechanics and heat transfer. These methods are 
(1) Experimental, (2) Theoretical, and (3) Numerical. The theoretical method 
is often referred to as an analytical approach while the terms numerical and 
computational are used interchangeably. 

In the experimental approach, a model is constructed and tested in a 
testing facility such as a wind tunnel. The flow variables, such as wall 
pressure and temperature can then be measured. The problem of producing 
required freestream conditions in the test section of this facility can be 
quite troublesome and time consuming. Since the facility, for example a wind 
tunnel, requires large amounts of energy for its operation, its operating 
costs are quite high. The experimental approach produces the most realistic 
answers for many flow problems; however, the costs are becoming greater 
everyday. 

In the theoretical approach, assumptions are made in order to simplify 
the problem. A closed form solution is generally sought. The main advantage 
of this approach is that general information which is usually in formula form 
can be obtained. Its disadvantage is that the problem is restricted to simple 
geometry and physics. 

In the numerical approach, a limited number of assumptions are made and a 
high-speed digital computer is used to solve the resulting governing fluid 
dynamics equations. The major advantage of this approach is that the problem 
is free of some of the constraints imposed on the experimental or theoretical 
approach. Thus, the numerical approach has the potential of providing 
information not available by other means. However, the approach does have 
some disadvantages. The storage and speed of present available computers pose 
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the limitations on the method. The other limitations arise due to the 
inability to understand and mathematically model certain complex phenomena. 
In spite of these limitations, the numerical approach is becoming more 
popular. The developments of supercomputers and the reduction in 
computational costs have made the approach appealing. 

A new methodology for attacking the complex problem in fluid mechanics 
and heat transfer has been developed and has become known as Computational 
Fluid Dynamics (CFD). Some of the ideas of this numerical approach are very 
old. Good surveys on the approach can be found in [1-13]. Also, the 
foundations upon which the whole field is built are now reasonably well 
covered in text books [14-24]. CFD is a science of producing numerical 
solutions to a system of partial differential equations which describe fluid 
flow. CFD is done by discrete methods and the purpose is to better understand 
qualitative and quantitative physical phenomena in the flow which then is 
often used to improve upon engineering design. CFD brings together a number 
of different traditional disciplines: fluid mechanics, the mathematical theory 
of partial differential equations, computational geometry, numerical analysis, 
and the computer science of programming algorithms and processing data 
structures. 

In the CFD calculation, the continuum problem of the differential 
equations is projected to some finite-dimensional space for the dependent and 
independent variables and then by solving the resulting discrete equations for 
the final set of numbers. Thus, the first step in CFD is to discretize the 
domain of the flow by laying out a network of points situated at a finite 
number of different locations of the independent variables, i.e., to create a 
grid. This brings about a so-called "grid generation" procedure. A "grid" is 
conventionally defined as a set of grid points in a coordinate system. The 
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word grid and mesh are also used interchangeably. Grid generation is an 
essential procedure in CFD. Accuracies and stabilities of the CFD 
calculations depend a great deal on the properties of these "grids". Grid 
points are generally generated by letting some coordinate lines coincident 
with the boundaries of the domain. The purpose of generating these so-called 
"boundary fitted" coordinates is to be able to apply boundary conditions 
directly when partial differential equations are solved on such grid. It is 
important that the boundary condition be represented accurately since the 
region in the intermediate vicinity of solid surfaces is generally dominant in 
determining the character of the flow. The procedure for generating a 
boundary-fitted coordinate can be divided into two catagories. They are 
partial differential equation methods and algebriac methods. In the partial 
differential equation methods, a set of partial differential equations, 
subjected to some boundary conditions, are solved to obtain a set of grid 
points. The partial differential equations may be elliptic, hyperbolic or 
parabolic. The partial differential equation methods offer the smoothness to 
the resulting grid points but generally require large amounts of computational 
time. The algebraic methods are based on mathematical interpolation functions 
and do not require the solution of differential equations or the use of 
complex variables. The primary advantages of algebraic methods are speed and 
directness. Regardless of how grid points are generated, all CFD calculations 
are usually done on a rectangular domain with a square grid. This is done by 
transforming the set of partial differential equations of interest, and the 
associated boundary conditions, to the curvilinear system. The grid points in 
the physical domain are, thus, mapped into a set of equally spaced grid points 
in a rectangular region called computational domain. With the transformation, 
the CFD calculations can be performed entirely on the fixed rectangular space. 
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regardless of the geometry or motion of the boundaries. Thus, numerical grid 
generation is the process of establishing an ordered and strategic distribu- 
tion of grid points in a physical coordinate system corresponding to a uniform 
distributions of grid points in a rectangular computational coordinate 
system. Detail on grid generation procedures are described in Chap. 2. 
Surveys of the methods including textbooks on grid generation procedure can 
also be found in references [25-30]. 

It is obvious that a grid which maps the entire physical domain onto a 
"slab" in the computational domain is very desirable. This type of grid, 
called single grid or single-block grid, offers the simplicity of the 
computational domain to the CFD calculation. However, for flow about complex 
configurations, the generation of a smooth and efficient single grid is very 
difficult. In some cases, especially those of configurations with several 
components, it may not be possible to obtain this type of grid at all. An 
alternative approach to this problem is to use several grids at once, each in 
a different coordinate system. The entire physical domain is, thus, 
subdivided into several subdomains. The generation of grids in different 
subdomain is generally independent from each other. This approach called 
"multiple grids" or "zonal grid" approach can be categorized into two groups: 
grid patching and grid embedding. For the patched grid approach, the 
subdomain grids are joined or patched together along common boundaries. The 
subdomain grids are overlapped rather than joined in the grid embedding 
approach. Multiple grids approach is becoming more common as the complexity 
of the configuration being considered in CFD is increased. However, the 
approach results in new boundaries which are not the physical boundaries. 
Even though, the solution procedure is done separately in each subdomain, 
certain boundary conditions are needed at these fictious boundaries. The 
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difficulty of the approach, thus, lies in the treatment of these boundary 
conditions. Since these boundaries are either joined or overlapped with other 
subdomain (or subdomains), some information needs to be transferred between 
subdomains so that the computation of the entire physical domain is 
consistent. The interpolation of flow variables between subdomains seems to 
be the simplest choice. However, this procedure does not result in a 

computational scheme which is conservative. The conservation of a 
computational scheme is important when the flow being considered contains 
discontinuities such as shock waves. A computational scheme is said to be 
conservative when it maintains the discretized version of the conservation law 
(conservation of mass, momentum and energy) exactly (except for round-off 
errors) for any grid size over an arbitrary finite region containing any 
number of grid points. For the multiple grids calculation the information 
needs to be transferred conservatively between subdomain grids. It has been 
suggested that fluxes rather than flow variables be transferred between 
subdomain grids, so that the resulting computational scheme is conservative. 
It can be shown that the conservation is easier to enforce in the grid 
patching approach than in the grid embedding approach. This study follows the 
grid patching approach along with the conservative treatment at the 
interfaces, i.e., places where two or more subdomain grids are joined 
together. The procedure is discussed in detail in Section 2.5. 

The viscous Navier-Stokes equations are the ultimate equation to be 
solved in CFD. However, Navier-Stokes flow simulation are presently still in 
the stage of research. A success in Navier-Stokes flow calculations relies 
not only on the numerical methods but also on the turbulent modeling. Since 
this study focuses on the concept of multiple grids approach, the inviscid 
Euler equations are sufficient to be used as the model equations. The Euler 
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equations result from dropping the viscous terms on the Navier-Stokes 
equations. Thus, the problems of computer storage and computational time 
including the uncertainties of turbulent modeling that arise in Navi er-Stokes 
calculation can be eliminated. Solutions of the Euler equations, though 
inviscid, give the correct phenomena to many of flow problems. The integral 
form of the Euler equations is applied to this study. The integral form may 
be important for the correct capturing of discontinuities in the flow. The 
discussion on the Euler equations is given in Sect. 3.1. The Euler equations 
are discretized by means of the centered finite-volume method. The finite- 
volume formulation is obtained by applying the integral form of the Euler 
equations to each grid cell of a given grid. The finite-volume method is a 
cell-oriented rather than grid points oriented. The main advantage of the 
method is that it can be applied to the general geometry without the need for 
a global coordinate transformation and it can tolerate the grid singularities 
since the flow equations are balanced only within the cells of the grid. The 
steady state solution is obtained by means of the time-dependent technique. 
The time derivative terms are reintroduced to the Euler equations and the 
steady state solutions are reached by explicitly marching the solution 
procedure in time from the initially guessed solutions. The three-stage 
Runge-Kutta integration scheme is used to serve this purpose. Since the 
transient solutions are of no concern, the local time step scaling is applied 
to accelerate the solutions to the steady-state. The linear and non-linear 
artificial dissipation terms are also added to the discreted Euler 
equations. The purpose of adding these terms is to impose an entropy 
condition which is required to eliminate the non-physical shocks. Also, the 
addition of the artificial dissipation terms helps to eliminate the 
oscillation of solutions which prevents the solutions from reaching the steady 
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state. Boundary conditions are of four types: solid wall conditions, 

inflow/outfl ow conditions, interface conditions and coordinate cuts. The 
finite-volume discretization along with numerical time integration and 
boundary conditions are described in detail in Chap. 3. The concept of local 
time step scaling and artificial dissipation are also addressed. 

The application of the approach to the flow over sphere at low Mach 
number and to the flow over a slender body at the supersonic Mach number are 
discussed in Sec. 4.1. Results obtained from these applications are very 
encouraging. Section 4.2 describes the application of the multiple grid 

approach to the flow over a Butler-Wing configuration. Again, good results 
have been obtained. The next step is to apply the approach to the 
internal/external flow on a fighter-aircraft. 
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2. GRID GENERATION 

2.1 Introduction 

Grid generation is an important procedure in CFD calculation. The word 
"grid" is generally used as a label for a complete set of grid points. Even 
though it is felt that grid generation and solution procedure are separate and 
distinct operations, in practice, these two operations can never be totally 
independent. This is because the accuracy of solutions depends upon grids on 
which the partial differential equations are solved. In turn, the logistic 
structure of the data (such as grid spacing), the location of outer boundaries 
and the nature of coordinate cuts are influenced by the nature of solutions. 
Perhaps the greatest problem of grid generations is not how to construct 
grids, rather, the problem is defining in sufficient detail what qualities and 
properties in a grid are desirable for a particular numerical method. 

The representation of boundaries is best accomplished when the boundary 
is such that it is coincident with some coordinate line, for then the boundary 
can be made to pass through the points of a grid constructed on the coordinate 
lines. Different expressions at, and adjacent to, the boundary may then be 
applied using only grid points and the intersection of coordinate lines, 
without the need for any interpolation between points of the grid. The 
avoidance of interpolation is particularly important for boundaries with 
strong curvature or slope discontinuities, both of which are common in 
physical applications. Likewise, interpolation between grid points not 
coincident with the boundaries is particularly inaccurate with differential 
systems that produce large gradients in the vicinity of the boundaries, and 
the character of the solution may be significantly altered in such cases. In 
most partial differential systems the boundary conditions are the dominant 


influence on the character of the solution, and the use of grid points not 
coincident with the boundaries, thus, places the most inaccurate difference 
representation in precisely the region of greatest sensitivity. The 

generation of a curvilinear coordinate system with coordinate lines coincident 
with all boundaries, so-called "boundary-fitted" coordinate system (Fig. 2.1), 
is thus an essential part of a general numerical solution of a partial 
differential system. 

Any partial system can be solved on the boundary-fitted coordinate system 
by transforming the set of partial differential equations of interest, and 
associated boundary conditions, to the curvilinear system. Since the 
bound ary- fitted coordinate system has coordinate lines coincident with the 
surface contours of all bodies presented, all boundary conditions can be 
expressed at grid points, and normal derivaties on the bodies can be 

represented using only finite difference between grid points on coordinate 
lines, without need of any interpolation, even though the coordinate system is 

not orthogonal at the boundary. The transformed equations can then be 

approximated using finite difference expressions and solved numerically in the 
transformed plane. Thus, regardless of the shape of the physical boundaries, 
and regardless of the spacing of the finite grid in physical field, all 
computations can be done on a rectangular field with a square mesh with no 
interpolation required on the boundaries. Moreover, the physical boundaries 
may even be time dependent without affecting the grid in the transformed 
region. Another major advantage of using boundary-fitted coordinates is that 
the computer software generated to approximate the solution of a given set of 
partial differential equation is completely independent of the physical 
geometry of the problem. Numerical grid generation is thus the process of 
establishing an ordered and strategic distribution of grid points in a 
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physical coordinate system cooresponding to a uniform distribution of grid 
points in a rectangular computational coordinate option. Some of the basic 
ideas of the use of boundary-fitted curvilinear coordinate systems in the 
numerical solution of partial differential equations are discussed in [31]. 
Figure 2.2 illustrates physical domain vs. computational domain. 

Two primary categories for arbitrary coordinate generation have been 
developed. They are algebraic methods and partial differential equation 
methods. The algebraic procedures include simple normalization of boundary 
curves, transfinite interpolation from boundary surfaces, the use of 
intermediate interpolating surfaces, and various other techniques. The 
partial differential system may be elliptic or hyperbolic. Included in 
elliptic systems are both the conformal, and quasi conformal mappings, the 
former being orthogonal. Orthogonal systems do not have to be conformal, and 
may be generated from hyperbolic systems as well as from elliptic systems. 

Algebraic transformations are attractive in that no numerical solution of 
a partial differential system is involved. Thus, the primary advantages of 
algebraic methods are speed and directness. Major disadvantages of these 
methods is the lack of smoothness that results when an elliptic partial 
differential system is used to generate the grid and truncation errors may be 
significant in regions where the grid is not smooth [32]. For instance, the 
results of Shang [33] show kinks in the solution corresponding to regions of 
rapid grid spacing change radiating outward from the boundary. It should be 
noted that local controls in the multisurface transformation [34] can be used 
to prevent nonsmooth boundary behavior (e.g., slope discontinuities) from 
propagating inward. Transfinite interpolation described by Gordon and Hall 
[35] in the early 1970 ' s is a highly generalized algebraic grid generation 
method. Transfinite interpolation is applied through a series of univariate 
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interpolations where blending functions and the associated parameters (point 
position and/or derivatives) determine a grid. For aerodynamic applications, 
Eriksson [36] and Rizzi and Eriksson [37] have adopted the original 
transfinite interpolation formulation to use only exterior boundary 
descriptions and derivatives at certain boundaries. They have also 
incorporated exponentials into the blending functions to concentrate the grid 
near an exterior boundary. The multisurface method [30,34] developed by Peter 
Eiseman provides formulas for grid definition based on grid descriptions of 
two boundary surfaces and an arbitrary number of intermediate control 
surfaces. Choosing interpolants (defined similar to blending functions) and 
the placement of the control surfaces determines grid shape and spacing. The 
multisurface method has been used by Eiseman in numerous applications [38,39] 
but most notably for computing grids about turbine cascades. The two-boundary 
technique [40-42] is based on the description of two exterior boundaries and 
the application of either linear or hermite cubic polynomial interpolation to 
compute the interior grid. For cubic interpolation, surface derivations 
combined with magnitude coefficients control the orthogonality of the grid at 
and near the boundaries. 

For the partial differential equation methods, a set of partial 
differential equations must be solved to obtain a coordinate system. The 
partial differential equations may be elliptic, hyperbolic or parabolic. The 
methods based on elliptic partial differential equations are more general 
(since all boundaries can be specified), and more fully developed. Typically, 
a pair of Laplace equations is solved subject to Cauchy-Riemann boundary 
conditions. The earliest successful development was formally reported by 
Winslow [43], who started with a Laplace system subjected to Dirichlet 
boundary conditions. Thompson, et al . [44] added periodic boundary conditions 
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to produce branch cuts for various topological configurations and who also 
suggested that control over the grid could be accomplished by altering the 
original Laplace system. The alteration is to consider a pair of Poisson 
equations by including specifications for the right-hand sides. These are 
called forcing terms and are general functions of the curvilinear variables. 
The particular form to be used was established later by Thompson et al . [45]. 
Without forcing terms, Mastin and Thompson [46] were able to show that the 
two-dimensional system analytically defined nonsingular transformation. 
Conformal mapping methods can also be included in the elliptic methods. Mehta 
and Lavan [47] have given a solution about a modified Joukowski airfoil 
accomplished by generating a coordinate system with a conformal Joukowski 
Transformation and solving the Navier- Stokes equations on the system. More 
examples of conformal mapping methods are given in Sampath [48], Wu et al . 
[49], and Napoli tano et al . [50]. When only one physical boundary is 
specified, hyperbolic partial differential equations may be used to obtain a 
grid by spatial marching from the given boundary. The remaining boundaries 
are determined by the solution and are geometrically unimportant in cases such 
as the external flow about a single object. A fundamental development has 
been given by Starius [51], and one which was well suited to body concavity 
has been presented by Stager and Sorenson [52]. The parabolic system can be 
applied to generate the grid between the two boundaries of a doubly-connected 
region with each of these boundaries specified [53-55]. The drawbacks of the 
hyperbolic scheme are: 1) the outer boundary can not be specified. 2) the 
scheme tends to propagate singularities of the boundary condition into the 
flow domain and 3) the solution may become unstable unless an artificial 
viscosity term is adequately added to the equations. On the other hand, the 
major dravA>ack of the parabolic scheme is that maintaining orthogonal ity of 
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grid needs much effort. Nakamura and Suzuki [56] have contained these two 
schemes into a single scheme that takes advantages of the two but eliminates 
the drawbacks of each. Both hyperbolic and parabolic methods have the 
advantage of being generally faster than elliptic methods, but are applicable 
only to certain configurations. 

It has been shown that the partial differential equation approach 
produces the smoothest grids for general boundary point distributions while 
the algebraic approach is the fastest procedure. Regardless of which approach 
is taken, creation of a computational grid requires (i) defining an accurate 
mathematical description of all solid surfaces in the computational domain and 
ii) generating an "appropriate" grid around these surfaces according to some 
criterion, usually with a specified point distribution. Graphic facility is 
very useful when three dimensions are involved. An important feature is the 
ability to rotate and translate grid surfaces in real time for inspection. 

In this study, algebraic approach has been taken due to its speed and 
directness. Two boundary technique [40-42] has been used to obtain boundary 
surfaces. The interior grid points have been obtained by applying the 
transfinite interpolation technique. Some detail of transfinite interpolation 
is given in the next section. Details of two boundary technique can be found 
in Ref. [40]. 


2.2 Transfinite Interpolation 

The idea of using interpol ation as a means of constructing general 
curvilinear coordinate systems stems from the fact that in most cases, the 
coordinates or grid points are known on several or on all of the boundaries of 
the computational domain and the problem consists of extending this grid into 
the interior of the domain. Interpolation from the boundaries into the 
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interior of this region can be accomplished by the so-called transfinite 
interpolation concept (sometimes referred to as the blending function 
method). Transfinite interpolation is a highly generalized algebraic grid 
generation method. Transfinite interpolation is applied through a series of 
univariate interpolations where blending functions and the associated 
parameters (point position and/or derivatives) determine a grid. The concept 
was originally developed by Coons [57] and subsequently extended by Gordon 
[58]. One of the earliest 2D grid generation applications using transfinite 
interpolation is described in Gordon and Hall [35]. A few examples of 3D 
applications are the works of Gerhard [59], Anderson and Spradley [60], and 
Spradley et al . [61]. In these applications, the transfinite interpolation in 
its simplest form is used, i.e. with no control of the normal derivatives of 
the grid coordinates at the boundaries. Eriksson [36,62] and Eriksson and 
Rizzi [63] have constructed a scheme which allows for the specification of any 
number of normal derivatives of the grid coordinates on the boundaries. The 
precise control of the resulting coordinate system or grid that this feature 
provides has made it possible to generate grids of advanced type that are both 
smooth and efficient in terms of resolution for a given number of grid 
points. 

Apart from giving good grid control, the transfinite interpolation 
concept offers speed and simplicity when implemented on computers. The speed 
factor is very important for 3D applications because the generation of a 
desirable grid for a given geometry is usually a process involving a series of 
grid generation runs with visual checks and adjustments in between. This fact 
is not always appreciated when evaluating the cost of grid generation. 
Naturally, a good graphics software package is an integral part of any 3D mesh 
generation system. 
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The theory of transfinite interpolation is a very general concept of 
mulitvariate interpolation and is outlined briefly. Let f(u,v,w) = [x(u,v,w), 
y(u,v,w), z(u,v,w)] denote a vector-valued function of three parameters u,v,w 
defined on the region u^cu <Up, v^cv <Vq, w^<w <w p . This function is known 
only on certain planes in the region, Fig. 2.3, 
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is needed to interpolate between these given planes. 

The transfinite interpolation procedure then gives the interpolated 

function f(u,v,w) by the recursive algorithm 
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The function f now defines a transformation from the region 
U|<u<Up.v^<v<Vq.Wj<w<w r in u,v,w space to some arbitrarily shaped region in 

■f 

the x,y,z space. It can be verified that if the specification of f on the 
planes u = u^,..«,u ; v = v ^*** v q and w = w^,.»», w p is continuous at the 
intersections of these planes, the explicit order of the interpolation 
directions chosen in the three-step algorithm does not affect the interpolant. 

The interpolation procedure just described can give any degree of control 
if a sufficient number of internal surfaces are specified, but the control is 
generally poor if no internal surface is defined at all. In order to improve 
the control while maintaining minimum input geometry data, a generalized 
transfinite interpolation procedure which uses derivatives of the 
function f in the out-of-surface direction, in addition to the function 
itself, can be defined. The effect of specifying out-of-surface derivatives 
of f (Fig. 2.4) is to introduce a direct control of the essential properties 
of the mapping function in the vicinity of the surface. The specified data 
are written as 
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which is simply the specification of f and on finite number of out-of-surface 

■f 

derivatives of f on the outer surfaces of the region u^u^u^, v^v^, w^<w<w^ 
in u,v,w space (Fig. 2.5). To interpolate f into the interior of this 
parametric box, a new set of univariate blending functions is defined as 
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The generalized transfinite interpolation algorithm is then written 
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The function f now defines a transformation from the region u 1 <u<u 2> 
v l <v<v 2’ w i <w<w 2 in UjV,w s P ace t0 some arbitrarily shaped region of x,y,z 


space. 
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Generally speaking, the method of transfinite interpolation, is a very 
simple and straightforward concept that offers virtually unlimited 
possibilities; but for any particular application, it is necessary to supply a 
certain amount of geometry data to obtain a certain degree of control of the 
transformation. It is up to the user to balance the requirements of minimum 
input geometry data and maximum control. 

2.3 Mapping Type and Singularities 

It is clear from previous sections that the first stage of grid generation 
procedure is the specification of grid coordinate data on the boundaries. 
Thus, the correspondence between the boundaries in the physical domain and the 
computational domain has to be made clear. It is then, necessary to determine 
the overall structure of the mapping between these domains. For a given 
geometry, there are generally several possible mapping types with different 
characteristics in terms of efficiency, coordinate cuts, singularities, etc. 
For example, there are at least six natural mapping types for the exterior 
region of a typical airfoil Fig. 2.6. All of these alternative mapping types 
give boundary-fitted coordinates but varies markedly in terms of grid 
efficiency, i.e. the resolution per grid point. It has been shown that the 
mapping type designated 0-0 is the most efficient for such a configuration 
[64]. The notation 0-0 is to be interpreted as "type 0 in the chordwise 
direction, type 0 in the spanwise direction", using the 2D notation shown in 
Fig. 2.6. Figure 2.7 illustrates the 0-0 mapping type of a wing-fuselage 
configuration. As shown in the figure, the entire wing is mapped to the bottom 
of the computational box, the entire outer boundary is mapped to the top and 
the combined plane of symmetry and fuselage is mapped to one of the side 
surfaces. The remaining three surfaces of the computational box constitute 
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cuts, i.e. they correspond to interior surfaces in the physical domain across 
which the various flow properties are continuous. 

Figure 2.8 shows that the 0-0 mapping type gives rise to two singular 
lines extending from the two tip corners of the wing to the outer boundary. 
Grid singularity is defined as a place in the physical domain where the 
Jacobian of transformation is zero or unbounded (depend upon how the Jacobian 
is defined). Grid singularities are undesirable but very often unavoidable and 
any practical finite difference scheme must be able to cope with them. If the 
physics near a singularity is not of primary interest, the finite difference 
solutions can be obtained in this region providing that the singular points 
themselves are excluded. Another practical approach to dealing with 

singularities is to leave the boundary surfaces open (Fig. 2.9). However, this 
requires that assumptions be made about the physics that must be included in 
the solution procedure. This study follows the so-called finite-volume method, 
which is a conservative cell-oriented method, and can be shown to be stable 
regardless of the type of singularity involved. The discussion on the finite 
volume approach is given in section 3.1. 

Since singularities always associate with the mapping types, and some 
types of singularities are more severe than the others, it is important to seek 
the best type of mapping for a given geometry, both from the viewpoint of efficiency 
and from the viewpoint of accuracy. For example, the C-H mapping has been the most 
popular type for flow computations around wings, even though this mapping has the more 
severe kind of singular line along the wing tip. Also, reference [64] has shown that 
this type of mapping is not as efficient as other types of mapping (for example 0-0 
type). This fact can be explained that the C-H mapping can be obtained by a simple 
"stacking" of 2D chordwise grids (c-types) in the spanwise direction, i.e. by a "quasi- 
3D" method. From the discussion in this section, it may seem that the price to be paid 
for using such a simple grid generation technique is high. 



20 


2.4 Multiple Grids 

The discussion so far have been limited to the topic of single-block 
grids, i.e., grids that map the physical domain onto a "slab" in the 
computational domain. This type of mapping is very desirable due to its 
simplicity. However, for very complicated geometries it can be difficult to 
generate single-block grid that are both reasonably smooth and efficient. An 
example of such a complex region is the exterior of a complete airplane with 
several lifting surfaces. Each component of an aircraft, in general, requires 
a grid system that is usually incompatible with the grid systems of the other 
components. Thus, the generation of a single boundary-fitted grid for the 
entire configuration is a difficult task, if it is possible at all. In such a 
global grid, control of grid point distribution, skewness and clustering will 
be difficult to achieve. For example, a grid which provides sufficient 
resolution of grid points in a region may result in an excessive number of grid 
points in other regions. Convergence to machine zero may not be achieved if 
the number of grid points is excessive. To simplify this problem, it is 
becoming more common to use several grids at once, each in a different 
coordinate system [65], for example see Fig. 2.10. This approach, called 
"multiple grids" or "zonal grids" approach (the terms "zone" or "block" is also 
used interchangeably), falls into two catagories: grid patching and grid 
embedding. The approach subdivides a complicated domain into subdomains which 
can accomodate easily generated grids. For the patched grid approach, the 
global grid is formed by patching together all the individual grids. The 
computed grid lines in adjacent grids may be made to align at the grid 
interface with complete continuity [66,67], or with continuous lines slope 
[64,68], or discontinuity in slope [69], or perhaps not align each other at 
all. Robert and Lee [69] combine the subdomain grid in such a manner that the 
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resultant global grid is continuous across patch boundaries. Moreover, grid 
irregularities frequently occur at the corners of the subdomain and at surface 
perimeter lines. Such irregularities impose constraints upon the choice of the 
numerical algorithm used for solutions of the flow equations. Lasinski et al . 
[70] have demonstrated a patched grid technique for solution of the thin layer 
Navier- Stokes equation. They solve the flow equations on each grid 
separately. The solutions are coupled by the transfer of boundary data at the 
coincident boundary points between grids. References [71-74] also illustrate 
the use of patched grid approach. The grid embedding approach does not need 
common boundary between grids, but rather, the various subdomain grids are only 
required to overlap to provide communication among grids for flow solvers. The 
development and analysis of solution procedures on grid embedding approach have 
been studied by Starius [75,76], Kreiss [77], and Mastin and McConnaughey 
[78]. The practical application of overlapping grids to the solution of 
problems in computational fluid dynamics has been demonstrated by Atta [65], 
Thompson [79], Steger and Buning [80], and Benek et al . [81]. Steger et al . 
[82,83] have applied the grid embedding technique to an airfoil/flap in 
incompressible flow [82] and in subsonic, compressible flow [83]. Atta and 
Vadyak [84] have obtained a potential solution for a wing/nacelle geometry. 
These studies have demonstrated that the technique can be applied to subsonic 
flow. However, for the transonic flight regime Benek et al . [83] have found 
that their single trial solution resulted in an ill-defined shock wave at the 
grid boundaries and exhibited poor convergence. The studies by Dougherty [85] 
indicate that for a different grid geometry and algorithm, these problems may 
not be too severe. Figure 2.11 illustrates an example of grid patching v.s. 
grid embedding. Early efforts to predict multiple-component configurations are 
based on the transonic small disturbance formulation [86-88]. Efforts to 
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predict the flow field about a complete aircraft configuration using a single 
grid approach have been made by Yu [89]. However, the requirement of exact 
boundary-fitted grids along certain boundary lines is relaxed. Thus, the exact 
implementation of boundary conditions is not obtained. 

The multiple grid approach has a number of advantages. First, the 
difficulty in generating three-dimensional grids for different types of complex 
configurations can be eliminated. Second, the approach allows different types 
of grid topologies to be implemented in each subdomain in order for grids to be 
mesh-efficient, i.e., more grid points near a solid body or shock and less grid 
points elsewhere. Since it is well known that skewness, rapid volume 
variation, and large cell aspect ratios degrade the convergence rate of an 
algorithm, it seems plausible that the enhanced grid point control afforded the 
multiple grids apporach will also result in improved algorithm performance. 
Third, it may also be computationally efficient to solve different equation 
sets in the various subdomain grids, such as viscous Navier- Stokes near the 
body and inviscid potential in the outer field. Chaderjian and Steger [90] 
have demonstrated this idea by solving the euler equations in one zone and the 
dual potential equations in the other for the transonic flow over a lifting 
ai rfoil . 

A common difficulty to the multiple grid approach is constructing a proper 
scheme for information exchange among the different subdomain grids. The 
information exchange has to be not only consistent with the governing 
equations, but should also lead to a stable efficient scheme. This so-called 
"interface conditions" is required to guarantee the convergence to a weak 
solution of the governing equations if the algorithm converges. The multiple 
grid approach results in new boundary within each subdomain grid, i.e., at the 
interfaces of various grids. Since these boundaries are not the physical 
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boundaries, it is important to treat grid points on the interfaces with care in 
order to transfer information from one grid to another accurately. The most 
obvious procedure is to interpolate the solutions in one grid to provide 
necessary boundary data for another. Since the classical interpolation 
formulas were not derived with conservation properties in mind, their use in 
finite-difference approximation on multiple grids would result in the loss of 
an exact conservation property. Eberhardt and Banganoff [91] have shown that 
shock waves crossing embedded grid boundaries can become ill defined and 
convergence is generally degraded when the interpolation procedure is used. 
They have also shown that the characteristic approach is superior but suggested 
that the use of conservative properties would be most desirable. For the 
existence of solutions to certain systems of partial differential equations, 
some conservation laws must be satisfied accurately. The nonlinear nature of 
the equations of motion permits solutions with discontinuities such as shock 
and slip surfaces. In order that such discontinuities assume the right 
strength and physical location, it is imperative that the scheme used for the 
calculation be conservative [92]. In a multiple grid calculation, it is 
important that the interforces are also treated in a conservative manner so 
that the discontinuities can move freely across the interfaces [93]. The need 
for conservative grid interfaces is also illustrated in ref. [83]. 

The question of conservation when switching between two different grids or 
numerical schemes has been considered by several authors. Warming and Beam 

[94] have derived transition operators for switching conservatively between 
MacCormack's method and a second order upwind scheme. Hessenius and Pulliam 

[95] have applied this transition operator approach to derive the so-called 
zonal interface conditions; this however, resulted in a significant loss of 
accuracy at the zonal interfaces. Rai [96] has developed conservative zonal 
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interface conditions for zonal grids which share a common grid line, and has 
provided accurate calculations demonstrating the shock capturing ability of the 
zonal grids with a discontinuity crossing zones. Combier et al . [97] have 
analyzed the zonal -boundary problem for a system of hyperbolic equations and 
used the compatibility equations to develop a zonal-boundary scheme. 
Reasonably good results were obtained for transonic channel flow. However, the 
use of the compatibility equations results in a scheme that is not conservative 
and, hence, unsuitable for problems in which flow discontinuity move from one 
grid to another. Rai et al . [98] have presented results obtained metric 
discontinuous grids; the integration scheme used is the Osher upwind scheme. 
Reference [82] provides the results obtained on overlapping grids in conjunc- 
tion with the stream function approach. 

In patched grid approach, conservation can be easily maintained at the 
patch interfaces. The extra computing time that is required to implement the 
zonal boundary condition is less than what is required for overlapping grids. 
This is because the necessary interpolations, that effect transfer between 
zones, are performed in a reduced number of spatial dimensions for patched 
grids. A problem in three dimensions only requires a two-dimensional 
interpolation procedure. This reduction in the number of dimensions in which 
the interpolation is performed does not occur for overlapping grids. On the 
other hand, overlapping grids provide more flexibility in generating grids 
because there are fewer constraints on the choice of outer boundaries for the 
different zones. Other disadvantages of grid embedding approach, beside that 
of interpolation, are: i) it is difficult to maintain global conservation ii) 
the accuracy and convergence speed of the calculation seems to depend on the 
degree of overlapping of the zones and the relative size of each zone, thus 
introducing a certain amount of undesirable empiricism in the formulation. 
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This study follows the grid patching approach in which the interfaces between 
subdomain grids in three dimensions are patched as plane interfaces. It can be 
shown that global conservation can be easily maintained for this type of 
interface. The study follows the method for transferring a conserved quantity 
from one generalized mesh to another which was first described by Dukowicz 
[99]. Ramshaw [100] has suggested a procedure in doing so which is similar to 
the method of Dukowicz, but is simpler and more direct. A computer program 
following the Ramshaw' s procedure has been written and tested with various 
types of grids and variables. The program has been working well for simple 
test cases. The objective of this study is to establish whether or not this 
technique is feasible for applications to realistic aerodynamic configurations. 

The grid generation of multiple grid does not in principle differ from the 
generation of single grid. The complete grid is computed by first dividing the 
entire domain into several subdomain grids and then "filling in" one subdomain 
grid after the other by transfinite interpolation. Eriksson [101] and Erkisson 
et al . [102] have obtained good solutions for the inviscid flow around an 
airplane by applying this concept. There, slope continuity (C 1 continuity) 
between subdomain grids is obtained by using osculatory interpolation, i.e., by 
using derivative information as well as grid point locations in the 
interpolation. The approach used in this study is different from that in Refs. 
101 and 102. Although the surface must be common between two subdomain grids, 
there is no restriction on grid slope or density across interfaces. This 
offers a great flexibility to the generation of each subdomain grid. Detail on 
the treatment of the conditions at the interface is given in sec. 3.3. 
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2.5 Brief Discussion on Conservative Rezoning Algorithm 

A method for transferring a conserved quantity from one generalized mesh 
to another, when the volumetric density of the quantity is assumed to be 

uniform within each grid cell of the original mesh, is described briefly 
below. This method was first described in Ref. [99]. Reference [100] 

suggested the procedure in doing so, which is similar in spirit but simpler and 
more direct. A computer program following this procedure has been written and 
is working well for example grids and a wide variety of choice of variables. 

By far the most common type of generalized mesh is the arbitrary 
quadrilateral mesh, which is convenient to work with because it has the same 
simple topological and logical structure as a square or rectangular mesh. The 
basic concept of the algorithm is simple. Consider Fig. 2.12, two grid 
surfaces are overlapping each other in some fashion. The conserved quantities 
Qo.jj of the original grid surface (Ao^- is the area of each surface mesh) is to 
be transferred to another grid surface in which An^- is the area of each 
surface mesh, Qn^- is denoted as the transferred quantity in each of these 

later surface mesh. Thus, Qn^- can be computed by 

Qn 1j = £ = E ( %U )(An ok t > < 2 ‘ 3 > 


Where An Qkjl is the portion of the area An^- which contains in the area 

Ao^ and the summation is up to the number of the original surface meshes 

contained in An.^ . And g Qkjl = Q ok / A olu = Volumetric density of Q ok£ , is 

assumed to be constant. The task now, is to find An . „ and the number of 

oka 

original surface meshes contained in An^-. The area of the polygon in 2D plane 
is given by [103]. 
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where the summation is over all the sides of p, and ej? is either +1 or -1 
according as p lies to the left or right, respectively, of side s. The 
endpoint coordinates (xj y*) and (x^ y^) are considered to be associated with 
the side s and not with the particular polygon. The overlapped areas are 
polygons p whose sides are segments of the old-mesh lines and the new-mesh 
lines. The number of sides of each type, and total number of sides, will be 
different for different overlapped areas. Each side is common to two 
overlapped areas, the one on the left (L) and the one on the right (R), and 
these overlapped areas may be considered to be associated with the side. The 
objective is to apportion a conserved quantity Q, whose volumetric density q is 
considered uniform within each cell of the old mesh, into the cells of the new 
mesh. It would be inefficient and difficult, to automate in a computer, to 
naively sweep over the overlap areas directly. Instead, Ramshaw [100] suggests 
to evaluate the same contributions by sweeping over the sides or segments s. 
The side or segment s is any side or segment of the polygon (overlapped 
area). The coordinate of the two end points of side s are denoted by 

(*i ^ 1 ^ ®nd ( x 2 y 2) ' 

If the side s is a segment of the old mesh then the quantity Q in the new 
mesh cell containing side s is therefore to be incremented by an amount 
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If the side s is a segment of the new mesh then the contribution to cell 

N 1 s s s s 

on the left is A s = 2 q o ^ x l y 2 ~ X 2 * 1^ while the contribution to the cell on 

N 

its right is just - a $ , where q Q is the volumetric density of the old mesh 

cell in which side s lies. 

0 N 

Adding a s and a s for each of new mesh cell yields the quantity Q contained 


in each of the new-mesh cell. 


28 


2.6 Applicetion to Two and Three Dimensional Hyperbolic Equations 


The first step toward the application of the technique to the CFD 
calculation is to apply the technique to solve some partial differential 
equations. The hyperbolic equations have been chosen not only because of their 
simplicity but also because of their hyperbolic nature which is similar to the 
equations of motion in supersonic flow. The hyperbolic equations can be 
written aSj 

in two dimensional space q t + aq x + bq y - 0 (2.6) 

and in three dimensional space q t + aq + bq + cq = 0 (2.7) 

where a, b and c have been treated as Constances. 

If the initial conditions are given as q = f(x,y) and q = f(x,y,z), the 
exact solutions can be found as 

q = f(x-at, y-bt) 
and 

q = f(x-at, y-bt, z-ct) 

for the two and three dimensional space, respectively. 

In two dimensional case, the equation has been solved on a two dimensional 
grid system which is changed into another grid system at some time. The 
procedure described in the previous section has been used to transfer the flux 
(in this case q itself is "flux") from one grid to another. The three 
dimensional equation is more suitable as the model equation of the equations of 
motion. The domain is divided into two subdomains which are independent from 
each other. Again the technique described previously is used to transfer flux 
(which is aq, bq or cq depending upon how the interface is oriented) across the 
interface. Results have been compared with the solutions from single grid 
calculations. Satisfactory results have been obtained. Details of this study 
can be found in Ref. [104]. 
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3. APPLICATION TO THE AERODYNAMIC CONFIGURATIONS 

The ultimate equations to be solved in CFD are the viscous Navier-Stokes 
equations. However, since solving these equations on modern day computers is 
still quite time consuming, they are often reduced to a simpler form. 
Solutions to these simpler equations, namely, stream function formulation [105] 
full potential equations [106-109] and Euler equations [110-115], have been 
obtained. The stream function formulation retains the generality contained in 
the full Euler equations. However, it is limited to two-dimensional or 

axisymmetric flows, and made difficult by the fact that the density in the 
transonic regime is a double-valued function of the unknown stream function. 
The full potential equation has been used as a standard model and has proved to 
be a helpful tool in the design of aircraft. As with the stream function, the 
potential equation can be solved by efficient relaxation techniques, and 
requires storage of only a single variable. Furthermore, it permits the 
solution of three-dimensional as well as two-dimensional flows. The primary 
disadvantages are the limitation to isentropic and irrotational flows. The 
isentropic assumption implies that shock waves captured in the transonic regime 
must be limited in Mach number to a value less than 1.3. The i rrotational ity 
conditions requires a uniform incoming flow in two-dimensional situations, and 
a free vortex condition three-dimensional flows. The potential equation will 
admit the existence of discontinuities in the flow field. However, these 

discontinuities are isentropic shocks, which do not represent true physical 
shock waves because they do not satisfy the Rankine-Hugoniot conditions. These 
shocks will be approximately of the proper strength and will exist in the 
proper portion in the flow field if the Mach number of the flow approaching the 
shock is less than or equal to 1.3. 

In this study, Euler equations are used as the model equations. Methods 
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based on the Euler model are useful tools in CFD since they offer more realism 
than potential methods and yet are simpler and more economic than methods based 
on the Navier-Stokes equations. A number of efficient and reliable numerical 
schemes have been developed for the Euler equations [110-115]. Even though 
viscous terms are neglected, solutions to Euler equations agree well with the 
experimental results. Shock waves captured in this model agree with the 
Rankine-Hugoniot relations regardless of their strength. And, more 
importantly, vortex sheets and vorticity can also be captured as weak and 
genuine solutions. The applications of numerical methods to solve the Euler 
equations range from the study of flow field around military aircraft and 
missiles where shock waves are strong, to more complex non-uniform shear flows 
past wings. Details of the Euler equations are given in section 3.1. 

The Euler solution procedure is based on a center finite-volume scheme 
with explicit Runge-Kutta time stepping [116]. This type of scheme was first 
used by Jameson et al . [117], but the present scheme differs significantly from 
this original scheme, mainly in the definition of the damping terms and the 
farfield boundary conditions. It has been extensively tested in both two and 
three space dimensions, for three different Euler models (the full equations, 
the constant-stagnation enthalpy model, and the artificial compressibility 
model for incompressible flow) and for both aerodynamics and turbomachinery 
application [118-121]. The finite-volume scheme is described in section 3.2. 

In most instances the solution to the first order steady state equations 
is desired. The steady state Euler equations change character depending upon 
the local Mach number. In a totally supersonic flow some very efficient 
methods exist for their solution. The method of characteristics and a simple 
marching procedure are two common approaches. In subsonic domain, however, no 
generally accepted method has yet been devised for solving this system. One 
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approach used for subsonic or transonic flows is to reintroduce the time 
derivative terms to the equations. The resultant set of equations is 
everywhere hyperbolic. A steady state solution can be obtained by marching 
from some initial guessed flow field through time until an axymptotic steady 
state is achieved. The initial conditions give rise to perturbation waves 
which move through the field as the solution progresses in time. The Euler 
equations have no inherent dissipation and, therefore, these waves must either 
be radiated from open boundaries or absorbed by the addition of artificial 
damping terms. The second difference and fourth difference damping terms are 
added to the Euler equation. The fourth difference terms are global and linear 
whereas the pressure-controlled second-difference terms are non-linear and are 
only activated around shocks. Boundary conditions are mainly of four types: 
solid wall conditions, interface conditions, inflow/outflow (farfield) 
conditions and coordinate cuts. Sections 3.3 and 3.4 describe the damping 
terms and numerical implementation of boundary conditions, respectively, in 
detail . 

Generally, to reach a steady state solution requires a large number of 
iteratives and a long computational time [122]. Since only steady state 
solutions are desired, and true time accuracy is of no concern, the concept of 
local time stepping is used to accelerate the convergence to steady state 
solution. This concept is introduced in Sec. 3.6. The explicit three-stage 
Runge-Kutta integration scheme is also addressed in Sec. 3.5. 

3.1 Governing Equations 

The Euler equations describing three-dimension, unsteady and compressible 
flows in conservation form can either be written in the differential form 
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where 


aQ + aF , aG + jH = 
at ax ay 32 


(3.1) 


p 

P u 

Q = pv 
pw 
E 


pu 

pv 

pw 

2 

pu + p 

puv , G 

puv 
2 , 

= pv + p 

pUW 

H = pvw 

pUW 

pVW 

pw + p 

(E+ p) u 

(E+ p)v 

(E+ p)w 


and 

p = density 

u = x - component of velocity 

v = y - component of velocity 

w = 2 - component of velocity 

E = total energy = internal energy + kinetic energy (in the absence of 

2 2 2 

Potential Energy) = pe + 0.5 (u + v + w ) 
or in the integral form 


d 

— / Q dxdydz = 6 (n »F + n .G + n .H) ds = 0 (3.2) 

dt J ft ft x y z 


where ft = arbitrary finite region 

n x = x component of normal vector at the boundary of the region 

riy = y component of normal vector at the boundary of the region 

n z = z component of normal vector at the boundary of the region 

the perfect gas equation of state is used to define the mean pressure P via the 
internal energy e: 


P = (y- 1) pe 


(3.3) 
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where Y = specific heat ratio = 

An assumption has been made, in writing Eqs. (3.1) and (3.2), that the fluid 
has undergone no external forces. It can be shown that the system of 
conservation laws given by Eqs. (3.1) or (3.2) is hyperbolic [17]. Thus, Eqs. 
(3.1) or (3.2) can be integrated in time in order to achieve a steady state 
solution (if such a solution exists). Equation (3.1) can be obtained by 
dividing Eq. (3.2) by n and then shrinking n to a point. This leads to the 
system of the differential conservation laws valid at that point if the partial 
derivatives are continuous there. The integral approach may be important for 
the correct capturing of discontinuities in the flow since it formally does not 
exclude discontinuities from the interior of n. This study follows the 
integral approach in which the difference equation are written directly from 
the integral system. Therefore, the method is a cell concept rather than a 
grid-point concept. The discussion on the method is given in Sec. 3.2. 

The nonlinear character of the Euler solutions generally permits solutions 
with discontinuities (shocks) where the differential Eq. (3.1) is no longer 
valid. The equivalence between Eqs. (3.1) and (3.2) is restored by allowing 
weak solutions to Eq. (3.1). However, both equations can give rise to 
nonphysical shocks unless an entropy condition is added. A "small" amount of 
artificial viscosity is added to the inviscid model for this purpose [96]. 
This artificial viscosity should also mimic the physical viscosity and create a 
primary vortex for flow past a highly swept wing at angle of attack. Although 
secondary vortices brought about by viscous effect, on the leeward side of the 
wing are not modeled, their effects on the primary vortices are small [123]. 
The Euler equations admit solutions with distributed vorticity but do not in 
principle contain any mechanism for generating vorticity. Any vorticity in the 
solution must be introduced either by boundary conditions or by shocks. Due to 
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the extra entropy condition shocks will lead to an increase of entropy and 

therefore also generate vorticity (Crocco's theorm) [124]. If the boundary 
conditions at the inflow boundary are such that vorticity is implied, this 
vorticity will naturally be convected into the domain and eventually be 
convected out at the outflow boundary. Furthermore, a solid boundary with 
sharp edge can also generate vorticity since attached flow around such an edge 
gives rise to shocks and thus also vorticity. In principle, this mechanism 
would act as an "automatic" Kutta condition [125] for the flow around an 

airfoil with a sharp trailing edge. However, some numerical evidence prevails 
that the combination of numerical errors and artificial viscosity will then 
produce vorticity and thus force the flow to separate at the edge. Section 3.3 
discusses the artificial viscosity model in more detail. 

For a finite domain it is necessary to construct suitable boundary 
conditions such that the desired steady state solution is obtained. The theory 
of absorbing conditions [126] is used in its simplest formulation. By 
linearizing the equations locally along the boundary and computing the 

characteristic variables along surface normals, it is possible to give the 
physically correct boundary information while maintaining good absorption of 
the transient error waves. The latter property is especially important for 
channel flows where stationary conditions are usually more difficult to obtain 
than for external flows. A more detailed description of these absorbing 

boundary conditions as well as other boundary conditions is given in Sec. 3.4. 

3.2 Spatial Finite-Volume Discretization 

Development of a method to solve the 3D Euler equations has been 
made (37,116,127-129]. It is a time-dependent finite-volume approach that uses 
multistage explicit time integration schemes together with centered space 
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differences. Significant features of this approach are integral conservation 
form, important for the correct capturing of shock waves and vortex sheets, its 
amenability to very general geometry without the need for or global coordinate 
transformation, and its toleration of grid singularities because the flow 
equations are balanced only within the cells of the grid [130], and not at the 
nodal point. It has been found that the time-dependent Euler equations permit 
solutions in which the flows separate from the leading edge of a sharp delta 
wing at angle of attack, without the implementation of a Kutta condition 
[125]. In contrast, separated flows are obtained by space marching methods 
only if a Kutta condition is enforced. 

The simplest way to derive the centered finite volume spatial 
discretization is to apply the integral formulation, Eq. 3.2, of the Euler 
equations to each mesh all of a given grid (see Fig. 3.1), i.e. 

d * -*■ " + * * 

— / Q dxdydz +6 (n -F + n *G + n »H) ds = 0 (3.4) 

dt si x y z 

i ,j .k an. 

i ,J s k 


where n- . . = volume element (i,j,k) 

1 j J . K 

By the mean-value theorem, Eq. (3.4) becomes 


dQ 

VOL - i ,j ,k + 
i ,j ,k dt 


6 F + 6 G + 6 H =0 
I i,j,k J i , j ,k K i , j ,k 


(3.5) 


where 


6 I F i,j,k " F i + V 2 .J s k " F i- V 2 >J >k 


5 jS,j,k ' ®i , j+ V 2 ,k ' ® i , j - V 2 , k 

6^H . = H . .. i.-H. .. 1 . 

K i ,j ,k i ,j ,k+ V 2 t • k— V 2 


(3.6) 
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and 


F i = SIX , -F . + SIY. , . * 6 . , 

i + V2 s j s k i +1 /2.j s k i + V2 s J » k 1+72 >J>k 1+72 >J,k 


+ SIY 

c 

+ SIZ 


i + V 2 J > k * i +1 /2 ,j,k 


F . = SIX . .F . + SIY , »G. , 

1 - V 2 5 J 5 k 1- V 2 >J J k 1” V 2 5 J > k 1- V2 » J 5 k 1- V2 »J > k 


+ SIZ 


i- V2 s J >k * i- V2 J,k 


( 3 . 7 ) 


G , = SIX , .F , + SIY , .G , 

i » J+ V 2 >k i s j + V 2 > k i j J + V 2 jk i ,j+ 72 .k i,j +1 /2 > k 

+ SIZ , .H , 

i sJ + 72 >k i >J + 72 > k 


SIX 


• F 


j <t — 01A 1 T 4 ^ SIY 4 *(3 *i 

i ,j- 72 ,k 1 s j- 72 s k 1 > J- 72 >k l.j-^.k i , j - V2 . k 

+ SIZ 


* 

•G 


.H 


i >J- V2 >k i ,J- V2 3 k 


H T = SIX , .F , + SIY , .G , 

i , J j k+ V2 i s j s k+ V2 i,j,k+y 2 i , j ,k+ V2 1 ,j,k+^ 


+ SIZ 


i j j ; k+ V2 i,J ! k + 1 /2 


H , = SIX . .F , + SIY . .G , 

i »j ik- 72 i ,j ,k- 72 i,j,k-y 2 i,J s k-V2 i ,J t k- y? 

+ SIZ , .H , 

i >J »k- 72 1 1 j ,k- V 2 


Since 
i.e. [?(})] 


>J >k 
V 2 J >k 


is 


located in the center of the cell but the flux function, 
must be expressed at its surface, some form of local 


interpolation of the neighboring discrete values Q must be devised. The 
simplest, and perhaps most natural, function is 
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5 i ,J + V 2 ,k> 

[f( ^) ] i,:,k+V 2 ° f( > 1 K^i,.j,k+V 2 ) ’ 


Expressions similar to eq. (3.8a) are obtained for £ and H. The average 
operator, u, is defined as 


U I * 1 + V 2 . J , k 1/2 ^i+l.j.k * •’l.j.k* 
“J*i,j + V 2 ,k =1/2( *1,j+l,k + 
VlJ,k+V 2 =1/2<, N|j,k*l + *1,J,k> ' 


An alternative is to compute the flux fuction separately for each of the two 
neighboring dependent variables and then average the two results, i.e. 


[W)] HV 2 ,j,.k + “I f( 5 i + V 2 ,j,k> 
Cf( 5 )] i,j + V 2 ,k ■ u / ( 5 i, j+ V 2 ,k ) 

C f ( <J) ] i ^ J kt - u/(f, ,j ,k+ 1/2 ) 


(3.8b) 


Similar expressions are obtained for (a and H. 

If the flux function were linear, alternatives (3.8a) and (3.8b) would be 
equivalent. For non-linear flux, only scheme (3.8b) provides the correct jump 
in Q across the shock. Thus, in this study, each term in Eq. (3.7) is defined 
as 
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F i + V 2 J ik " 1/2 ^i+l.j.k + F i ,j ,k^ 

F i , j+ V 2 " 1/2 ^ F i,j+l.k + F i ,j ,k^ 

F i,j..k+ 1/2 = 1/2 (F i,J,k+l + F i,j,k ) ' 

Similar expressions are obtained for $ and ft. Finally, the other terms are 
expressed as 


SIX i+V 2 J,k 

V 2 ((y i+ 

V 2 a j- V 2 .k+ V 2 " 

*i + 

V 2 > J + V 2 »k- I /2 ^ 
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Z H 
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V 2 ,j- V 2 >k+ V 2 ” 
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SIZ i+V 2 ,j,k = 
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V 2 J- V 2 jk+ V 2 " 
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V 2 J+ V 2 »k- 1/2 ^ 


•(y i+ 

V 2 »j+ V 2 ,k+ V 2 " 


V 2 ,j- l/ 2 ;k- l/ 2 } 


-(y i+ 

V 2 J- V 2 ,k+ 1/2 ■ 

yn 

V 2 »j+ V 2 ,k- V 2 ^ 


* (x i + V 2 ,j+ V 2 ,k+ V 2 ' X i+ V 2 , J- V 2 ,k- V 2 } ^ 



SIX i ,j+ V2 ,k~ ^2 (( y i+ V2 > J+ V2 »k- V2 " y i- V2 ,j+ V2 s k+ V2 ^ 

*^ z i+ V2 ,j+ V2 .k+ V2 ’ Z i- V2 J+ V2 »k- V2 ^ 

~^ z i+ V2 ,j+ V2 ,k- V2 “ Z i- V2 ,J+ V2 ,k+ V 2 ^ 


^ y i+ V2 ,j+ V2 ,k+ V2 ” y i- V2 >J+ V2 ,k- V2 ^ ^ 

SIY i > J+ V2 ,k~ ^ ^ z i+ V2 , j+ V2 jk- V2 " Z i~ V2 > J + V2 >k+ V2 ^ 

*^ x i+ I/2 jj + V2 >k+ V2 X i- V2 >j+ V2 »k- V2 ^ 

^ x i+ V2 j J + V2 >k- V2 X i- V2 1 J + V2 >k+ V2 ^ 

• (z i+V 2 J+V2 jk+ V2 " Z i- V2 J+V2 ,k- V2 } ^ 

SIZ i > j + V2 ,h =1/ 2 ( (x i+V 2 J+V2 1 k - V2 " X i-V2 , j+ V 2 ) k+V 2 ) 

*^ y i+ V2 :J+ V2 »k+ V2 " y i- V2 > 0 + V2 »k- V2 ^ 

~^ y i+ V2 >j+ V2 >k- V2 ” y i- V2 ! J + V2 >k+ 1/2 ^ 

* (x i+ V 2 ,j+ V2 >k+ V 2 " X i- V 2 , 0 + ] /2 »k- V 2 } ) 

SIX i , J ,k+ 1/2 = 1/2 « y i- V 2 ,J+ 1/2 .k+ 1/ 2 “ y i+ I /2 ,j- 1/2 ,k+ I /2 } 

*^ z i+ V2 ,J + V2 »k+ I/2 ' z i- I/2 , j- 1/2 s k+ I/2 ^ 

”^ z i- V 2 ,j+ I/2 .k+ I/2 " z i+ V 2 jj- I/2 5 k+ 1/2 ^ 

*^ y i+ I/2 ,j+ V? ,k+ l / 2 " y i- I/? ,j- I/2 »k+ I/2 ^ ) 
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; j ; k+ V2 ^ ^ Z i- V2 >J + V2 s k+ V2 " 2 i+ V2 J- V2 s k+ V2 ^ 

*^ x i + V2 )J + V2 s k+ V2 X i- V2 jj- V2 s k+ V2 ^ 

^ X i - V2 > j+ V2 jk + V2 X i + V2 s J- V2 jk+ V2 ^ 

*^ z i+ V2 J+ V2 »k+ V2 ~ Z i- V2 ?j- V2 >k+ V2 ^ ^ 

^ IZ i ,j : k+ V2 ^ ^ x i- V2 J+ V2 ,k+ V2 ” X i+ V2 >J- V2 ,k+ V2 ^ 

*^i+ l / 2 ,j+ V2 >k+ V2 " y i- V2 J- V2 >k+ V2 ^ 

~^i- V2 ,j+ V2 ,k+ V2 " y i+ V2 >J- V2 ,k+ V2 ^ 

* (X 1+ V2 .J +1 /2 , k+ 1/2 " X i- V2 s J“ V2 ,k+V 2 }) 

From these formulas it is clear that the only quantities needed from the 
coordinate transformation are the x,y,z - coordinates of the grid points. 
Equation ( 3 . 5 ) together with ( 3 . 8 b) leads to a spatial-difference operator 
completely centered in all three coordinate directions, which is second-order- 
accurate in space if the variation in grid size is reasonably smooth. 

The finite-volume discretization bears some similarity to both the 
conventional finite-difference and finite element discretizations. Its 
formulation, like the finite-element procedure, begins with the integral 
equation. Its difference stencil is that of a finite-difference scheme, but it 
differs in that cell-averaged instead of point quantities are differenced, and 
this gives a significant distinction near a grid singularity. In the finite- 
volume formulation, the flux quantities can be defined and remain finite even 
in the presence of the grid singularity, since Eq. ( 3 . 5 ) is balanced in the 
interior of the cell where no coordinates are used. The usual grid-point 
methods may not yield this feature without special programming consideration. 
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Eriksson [130] has concluded that without any modification the finite-volume 
technique remains stable in the presence of a grid singularity, but its 
accuracy decreases to somewhere between first and second order in space. 
Without alteration the finite-difference scheme is unstable even if the 
singularities is straddled. However, if a limiting form of the difference 
scheme is derived at the singularity point and implemented in the computer 
code, stability of the finite difference scheme can be restored. The important 
aspect of the finite-volume approach is that it suits well with the 
conservative rezoning approach used in this study. This is true because fluxes 
are obtained as the average values at the center of the cell faces, i.e., no 
interpolation from grid points is needed. So, the order of the interpolation 
scheme does not play a role in the interface treatment. 

3.3 Artificial -Viscosity Model 

The central difference schemes to solve the Euler equations are inherently 
dispersive and not dissipative. Even for linear problems, central-difference 
schemes admit as a solution so-called sawtooth waves. The non-linear nature of 
the Euler equations gives rise to an aliasing phenomenon whereby these short 
waves interact with each other, vanish, and reappear as distorted long waves. 
In nonlinear transport there is a mechanism by which energy migrates from long 
wavelength motion to progressively shorter and shorter scales until it is 
removed from the flow by molecular viscosity. The Euler equations possess no 
such viscosity so, in the discrete representation, this energy would migrate to 
the smallest scale resolvable on the grid and then returns to large-scale 
motion via aliasing, which is non-physical and would make a steady state 
unattainable [131]. In general, these defects could be dealt by digital 
filtering techniques. However, further deficiencies arise. The nonlinear 
conservation equations admit non-unique weak solutions when shocks are to be 
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captured. An entropy condition has to be supplied in order to obtain the phys 
cally correct weak solution [132]. A standard way to invoke an entropy 
condition is to model the true physical process inside a shock by the addition 
of a small dissipation term to the convective differences. This so-called 
artificial viscosity mimics the real physical viscosity not only by involving 
an entropy condition but also by removing the short-wave motion out of the 
flow. 


A number of literatures has been developed on construction of such 

artificial viscosity models, but they vary in detail from method to method. 

The construction of the models is arbitrary except for classification according 

to its order of magnitude in terms of grid spacing. In this study, the 

dissipation is introduced at the same time as the transport process. Its 

magnitude lies in or below the range of the trunction error of the discrete 

approximation. The total difference operator ?(qT therefore consists of: (1) 

the convective part ? ($) that results from discretizing the Euler equations 

in space by the centered finite-volume scheme, and ( i i ) the dissipative part 

? n ($) • Thus, equation (3.5), can be written as 
D« 


d > 

— Q 
dt ijk 


F (Q ) 
C i jk 


+ F (Q ) = 
D ijk 


F(Q ) 

ijk 


(3.10) 


The total discrete dissipative operator F D (Q 1 -j k ) includes its own 

artificial boundary conditions, and comprises both linear and nonlinear terms 

according to ?($..,) = f(C..,) + D$. , where D is a constant matrix. The 

a D v i jk' ijk' ijk’ 

nonlinear expression ^ ( C j ^ ) is designed to provide dissipation at 
discontinuities, whereas the linear one is formulated to suppress spurious 
solutions (sawtooth waves) and to control the migration of energy from large to 
subgrid scales. 
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3.3.1 Nonlinear Artificial Viscosity 


The nonlinear artificial viscosity in the interior of the domain is 
expressed by 


f 

ijk 


■ *i s i ts i ( VV 


s.Cs (Q. ..)«,] + 

J J ijk J 


WV 


0)0..,, 

k ijk 


(3.11) 


where x is a constant and Sj.Sj and s k are coefficients that depend on the 
solution field through the pressure according to 


s i = K LP ijkf s j = 'hj LP ijk'l and \ = l J k LP ijk't' 


2 2 2 

where fit . 6j and 6 k are central difference operators, 


6 1 ^i j k ~ ^i+1 , j ,k " 2 ^i , j ,k + ^i - 1 , j , k 


6 J ^ijk “ *i,j+l,k " 2 ^i , j , k + ^i , j - 1 s k 

5 k ^ijk " *i,j,k+l " 2 ^i , j , k + ^i,j,k-l 

and LP ijk = Log (P iJjk ). 

These coefficients are normalized by their maximum value so that their 
magnitudes lie between 0 and 1. Their purpose is to sense non-smooth flow and 
increase the filtering of large gradients so that in effect an entropy 
condition is enforced. At the boundaries, the coefficients Sj , Sj and s^ are 
set to zero. 


3.3.2 Linear Artificial Viscosity Team 

At all interior cells, the fourth-difference operator is used and the 
linear artificial viscosity is expressed as 
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where 


DQ 1Jk - - y«4 * «J + «k> Q ijk’ 

s I*ijk ' *1-2, J,k • 4 *1-l,j,k + 6 * 

5 0*1jk = *1,j-2,k - 4 *i,j-l,k + 6 * 

S k*1jk = *1,j,k-2 - 4 *i,j,k-l + 6 * 


i ,j,k 
i ,j,k 

1 » J .k 


(3.12) 

4 *1+l,j t k + ^i+2,j ,k 
4 *i ,j+l,k + *1 ,j+2,k 
4 *1,j t k+l + , j ,k+2 5 


and y is a constant. Linear extrapolation is used at the boundary cells. 

3.4 Boundary Conditions 

For the computation of many fluid dynamic problems more difficulty is 
encountered in satisfying the boundary conditions than in balancing the 

differential equations at the interior points of the flow field. This is so 
because on the boundary not all of the flow variables are specified by the 

boundary conditions, and there remain more unknowns than equations. While 
transformation to a boundary-fitted coordinate system does reduce to one the 
number of unspecified boundary variables necessary for differencing the 

interior field, namely the pressure [133], still some way is needed to couple 
these unknown values of pressure to those in the interior in a manner 
consistent with the boundary conditions. Improper treatment of the boundary 
conditions can lead to serious errors and perhaps instability. In order to 

treat the flow exterior to a domain an artificial outer boundary must be 

introduced to produce a bounded domain. This is an artificial boundary in the 
sense that the actual flow in the physical domain is open, whereas, the 
computational space must, for practical reasons, be closed. The numerical 
conditions, therefore, should allow phenomenon generated in the computational 
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domain to pass through the boundary without undergoing significant distortion 
and without influencing the interior solution. Thus, the maximum amount of 
transient energy can escape from the field so that the time-dependent solution 
converges to the steady state. Engquist and Majda [124] have presented a 
mathematical theory for the practical application of local absorbing boundary 
conditions at artificial boundaries. Their 'First Approximation' is adapted to 
this study. 

Four distinct types of boundary conditions: conditions at solid walls, 
periodic conditions across coordinate cuts, flow into or out of the artificial 
boundary, and conditions at the interfaces, are discussed below. 


3.4.1 Solid Walls 


At a solid wall the mass flux is zero but the surface pressure contributes 
to the momentum flux, Eg. (3.4), thus is simplified to 


/ (n «F + n »G + n *H) ds = f Sds 
solid x y z J sol id 

wall wall 


(3.13) 


0 

A 

n »p 
x 

where S’ n -p 

y 

A 

n »p 
z 

0 

Equation (3.13) is used to derive the contribution from those mesh walls which 
conincide with a solid wall. For example, if j = V2 i s denoted for these mesh 
cell walls (Fig. 3.2) then Eq. (3.13) is approximated by 
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A -»• A •+ A -*■ 
f (n «F + n *6 + n »H) ds = 
J sol id x y z 

wal 1 


0 

SJX , »p , 

i, \ ,k i, V 2 ,k 

SJY , .p , 

i , V2 ,k i , V2 s k 

i , V2 ,k H i , V2 ,k 

0 


SJX i . 1/2 ,k = V 2 « y i + 1/2 , 1/2 ,k- 1/2 ' y i- 1/2 . 1/2 .k+ V 2 } 

* (z i+ 1/2 , 1/2 »k+ 1/2 " Z i- 1/2 . 1/2 »k- V 2 5 

■ (2 i+ 1/2 , 1/2 ) k- 1/2 ' Z i- 1/2 , 1/2 >k+ l/ 2 } 

*^ y i+ 1/2 > I /2 »k+ 1/2 " y i- I /2 , 1/2 t k- I /2 ^ ^ 


SJY i , V2 ,k = 1/2 ( (z i + l / 2 , 1/2 ,k- V 2 " z i~ V 2 , V 2 ,k+ V 2 } 

*^ x i+ I/2 , 1/2 ,k+ I/2 x i- I/2 , 1/2 ,k- I/2 ^ 

' (x i+ l / 2 , l / 2 »k- 1/2 ' x i- 1/2 , 1/2 »k+ I/2 } 

• (z i+ 1/2 , 1/2 >k+ 1/2 ' z i- 1/2 , 1/2 >k- I/2 } ) 

SJZ i , 1/2 jk = 1/2 ( (x i+ I/2 , l / 2 .k- I/2 " x i- l / 2 , 1/2 ,k+ V 2 5 

*( y i+ I/2 , l / 2 ,k+ I/2 " y i- I/2 , 1/2 ,k- I/2 ^ 

"^ y i+ 1/2 5 1/2 »k- I/2 " y i- I/2 , 1/2 , k+ V 2 ^ 

*^ x i+ I/2 , I/2 ,k+ I/2 " x i- I/2 , 1/2 ,k- I/2 ^ ) 


and P-j I/2 k a PP rox1mat ed by a linear extrapolation from the 
interior cells, i .e. . 


( 3 . 14 ) 


center of the 



3.4.2 Coordinate Cuts 


At a coordinate cut the physical space folds on to itself and the 
condition on the flow at the computational boundary is periodicity. This 

boundary, in fact, does not exist as a physical boundary. It is a boundary 

only in a practical programming sense. Thus, it does not influence the 

solutions in the interior. An example of this type of boundary can be seen in 
Fig. 3.3. 

3.4.3 Inflow/Outflow Boundaries 

As mentioned previously, these boundaries are artificial boundaries for 
practical reasons. The theory of absorbing boundary conditions [123] is 

applied to convert the transient energy out of the flow field so that the 

steady state solutions can be reached. This is done by linearizing the 
governing equation locally and computing the characteristic variables in the 
normal direction. Those characteristic variables which are advected into the 
domain are then fixed to the desired values whereas those which are advected 
out of the domain are linearly extrapolated from the interior to the 
boundary. The resulting complete set of characteristic variables is then 

transformed back to the primitive variables and used to compute the desired 

fluxes. The concept of Engquist and Majda's 'First approximation* is described 
below. 

A corresponding one-dimensional system of linear hyperbolic equation can 
be written in the characteristic form as 
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where $ represents the characteristic variables and the Jacobian matrix A can 
be written as 

0 age 

2 2 

ac - uU + kaV ct(l-k)u + U 0U - kaV eU - kaW 

2 2 

A = gc - vU + kgV a v - kgu g(l-k)v + U ev - kgw , 

2 2 

ec - wU + keV aw - keU gw - keV e(l-k)w + U 

■+ ■> 

where a = s • e ; e = s - e ; e = s • e 
x y ^ 

2 > -*• y-1 

U = aU+gv + eW ; V = V*V, K = 

Y 

c is the local speed of sound and s is the surface area of the cell face 
coincide with the artificial boundaries. 

The eigenvalues x of A can be found by solving det (A-xI) = 0, as 
Xj = U, x 2 = U, x 3 = U - a + . x 4 = U - a_ , 

1, ... r l . 2, ,2 2.2 2 2... \ 

where a + = V 2 kU + [- k U + c (a +g +e )] 

The left and right eigenvalues associated with these four eigenvalues make up 
the row and columns of the transformation matrices T ^ and T respectively which 
diagonalize Eq. (3.15): 



— $ + A — ^ = 0 A. 

8t 3X \o 



where 

c 

; = T _1 $ , A = T _1 AT = x 3 

U A 4 

• 

(3.16) 

After the 

intermediate variables 




U = 0U + aV , V = - eV + gw , W = eU - aW , 
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5 = a + 8 2 + e 2 , Q + = ku - a + , R + = e Q + - ?kw . p + = kwa + + e c 2 
have been defined for the sake of simplicication, it can be found that 


T = 


and 


/V 

kU 

kvll - g(kV 2 + c 2 ) 

- kuU + a(kV 2 + c 2 ) 
0 


0 

R 


R 

V 

+ 

uR + + a p + 


uR + ap 

IV 

w 

vR + + gp + 


vR + gp 

u 

wR + - k(aU + 

Bv)a + 

wR - k(au + gv)a 


- (a 2 + 

„ 2 )c 2 

- (a 2 + g 2 )c 2 


C V 2 - U 2 -eW + gU 

d l d l 

[kw(U 2 - 5V 2 ) [kw( - eU + e W) 
+ C 2 (eU - cw)]/d 2 -aec 2 ]/d 2 


eV - oU aW - gV 



[kw( a U - e V) [ — kw( aW-gV) 

-geC 2 ]/d 2 +(a 2 + g 2 )c 2 ]/d 




k?v -(U+a + )Q + 


-k£u + oQ+ 



-kgv + gQ+ 




k5v 2 -(U+a_)Q- 



-k£U + qQ- 


d 


4 


-kgv + gQ- 


d 


4 


The factor dj, d2- d0 and d 4 in the dominators are normalizing coefficients so 
that T’ 1 T equals the unit matrix. 

For the one-dimensional case it is well known that the number of 
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conditions to be imposed in a cell at the outer boundary should equal the 

number of characteristic directions that enter the computational domain. Four 

typical cases are shown in Fig. 3.4. With subsonic inflow the implementation 

is to set the three ingoing characteristic variables , <j>^ and <j>^ to 

(41 

their free-stream values, linearly extrapolate the fourth <j> v ' from the compu- 
tational field, and then solve for the original unknowns Q = T<j> . At outflow 

it is 4 /^ that is given the values of undisturbed flow, and 4 >^, 4 >^ 

and 4 >^ are extrapolated from the computational field. 

3.3.4 Interface Conditions 

An interface, here, is referred to as a common boundary where two or more 
subdomain grids are patched together. As mentioned previously, these 
boundaries are not physical boundaries and may compose of grids of different 
topologies. Care must be taken in order to treat these boundary conditions 
which exist because the computation is done on each subdomain grid 

individually. This study follows a conservative approach which offers the 
conservation of fluxes at the interfaces. The conservative treatment at the 
interfaces may be important for the correct capturing of discontinuities 
crossing them. The discussion on the conservative rezoning algorithm is given 
previously. An example of its application is described below. 

Consider the common boundary (interface) between two grids as shown in 
figure 3.5. The application of the finite volume approach to the cell 
(i,j, NK-1) requires the integrated fluxes, = hjj^ . A(j^ at the cell face 
k = NK - V2 » where a(^ is the cell surface area at k = NK - V2 • For the 
interior cell, fluxes at the cell walls are computed by taking the average of 
the flux functions evaluated at cell centers, i.e., h S!j,k + l/ 2 1 / 2 (h i,j,k 
+ h.. j k+1 ). At the cell face k = NK -V 2 0 f grid 1 ,h^ j which contain in 
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grid 2 are needed. They are computed by interpolating the flux functions 
evaluates at cell centers of grid 2. 

For grid 2, the integrated fluxes h(^. are required for the 

computations of cell k=l. They are obtained by applying the conservative 

rezoning algorithm described previously. To illustrate how the algorithm is 

applied, consider an interface between two grids as shown in figure 3.6. Here 

the subscript k is dropped since the interface lies on the cell wall of 

constant k. The conservative treatment at the interface requires that the 

(21 

integrated fluxes going through an area A u be the same for both grids, i. e., 
H kj^ = H k^ ' frroni Fl S • 3,6) H k^ can be evaluated as 




P 

E 

n=l 


H (1 >. 

1J 


a (2) 

ktij 

A (1) 

1J 


P 

E 

n=l 


M) *( 2 ) 

h . . • A. 
ij kAij 


(3.17) 


where A^| in the portion of the area A^ which contain in the area kiV and 
p is the number of the surface mesh of grid 1 which contain in the area 
a£ 2' . in this manner the conservative rezoning algorithm can be applied to 
compute 


H 


(2). H (l). ? h U). A (2) 

H k* J =1 h ij A ki1j 


k£ 


(3.18) 


since h(^ are known from the previous interface treatment for grid 1. 

’ J 
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3.5 Numerical Time Integration 

The complete semi -discrete scheme. Eq. (3.10), including all boundary 
conditions defines a unique system of nonlinear ordinary differential 
equations 

u t = F(u) ; u ( 0) = u Q (3.19) 

which must be integrated in time by numerical means. The explicit three-stage 
Runge-Kutta scheme presented by Gary [134] is used in this study. This scheme 
is defined by the following algorithm: 

u*(t„ + i) = u(t n ) + AtF(u(tn) ) 

u **(t n+ i ) = u(t n ) + V 2 AtF(u (t n )) + 1 / 2 4tF(u*(t r|+1 )) (3.20) 

u ( t n +l) = “(t n ) +1 / 2 StF(u(t n )) + V 2 4tF(u**(t n+1 )) 

It can be shown that this scheme, when applied to the semi -discretized Euler 
Eqs. (3.10) is second-order accurate and is stable with a CFL-number of at 
most 2. 

The multistage two level schemes of the Runge-Kutta type have the 
advantage that they do not require any special starting procedure, in contrast 
to leap frog and Adam Bashforth methods, for example. The extra stages can be 
used either (1) to improve accuracy, or (2) to extend the stability, region. 
Another advantage of this approach is that the properties of these schemes 
have been widely investigated, and are readily available in textbooks on 
ordinary differential equations. 

3.6 Local Time Step Scaling 

As mentioned previously, to reach a steady state solution by explicit 
methods, generally, requires a large number of iteration and a long 
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computational time. This is so because time step used in explicit methods is 
restricted to a maximum value according to the CFL (Courant-Friedrichs-lewy) 
condition [135]. The maximum time step is usually determined by the smallest 
grid spacing. On a highly stretched grid, this maximum time step can be 
extremely small. In the applications where only steady state solutions are 
desired, and true time accuracy is of no concern, it has been found [136] that 
the use of the "local time step" scaling is highly beneficial to accelerate 
the convergence to steady state solutions. The simplest way to derive this 
scaling is by a local Fourier analysis. This concept is outlined below. 

To obtain a better understanding of the concept, let's consider the Euler 
equations in two space dimensions: 

Q + F(Q) + G(Q) = 0 (3.21) 

t x y 

Assuming that the mapping x( 5 ,n), yU-n) between the physical x-y space and 

the computational £,n - space is smooth. Equation (3.21) can be transformed 

into the computational space as 

(JQ) + (y F - x G) + (- y F + x G) =0 (3.22) 

t n n 5 £ £ n 

where J = Jacobian of transformation = x y - x y 

£ n ty £ 

Equation (3.22) is integrated over the region D. . and can be written as 

1 J 

~ J/ n QJdcdn - 4 [(y F - x G)d 5 + (y F - x G)d n ] = 0 (3.23) 

dt D, 3 D 5 5 n n 

ij ij 

The computational space has been discretized according to 

£ i = £ 0 + i AC ; hj = o 0 + JAn 

and same notation have been defined as 

x. . = xU., n.) .. y. . = y(£., n.) , Q. . = Q(x. y. t) 

i »J i J i.J i J i,J i,J 

+ + + + + + 

F. . » F(Q. .) , G. . = G(Q. .) 

1 5 J i ,J 
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5 g^i,j ^i+ V2 : j V2 J ’ .J+ V2 J- V2 

U Q = V 2 (Q . 1 + Q. i, .) , w Q. . = V2 (Q . .1 + 5 . . 1. ) 

g i,j 1+ V2 jJ i - V2 s J n i s J i,J+V 2 i,J-V 2 


Equation (3.23) can be approximated by 


d t 
dt 


Q .// Jdgdn + 6 [- (6 y. .)(y F. .) + (5 x ,)(u G. .)] 
i ,j 1 u n L g i,j n g i,J n i .J 

i , j 

+ 6 [(6 y. .)(y F. .) - (5 X. ,)(y G. .)] = 0 
g n i ,j g i ,J n 1 5 J i >J 


(3.24) 


Equation (3.24) is a semi -di screte centered scheme for the two dimensional 
Euler equation. To derive the local time step scaling for equation 3.24, 
first the transformation matrics are freezed at mesh point (i,j), i.e. 


x (g,n) - x^(g.,n j .) = x ? ; x^(g, n ) - x^g.,^.) = x^ 

y 5 (g,n) - y^U, ,nj) ■ y K ; y n (c.n) - y T1 (c i s nj) = y n 
J(g.n) = x ? (g,n) y^g.Tj) - x^(g,n)y (g,n) - - x^y ? = J (3.25) 

Thus, Eq. (3.23) becomes 


J AgAn— Q. . - ! / 2 Agy (F. - F. ) + \ Agx (G . 

dt i,j g i ,J+1 ijJ-1 5 1 >J 


- G ) 
j+1 1J-1 


■ y ^ n J „ ( Vj - F l-l,j } - 1/24n \ (G m,j - Vl.J 1 = 0 

(3.26) 

Next, the flux functions F(g), G(q) are linearized around Q-jj) according to 


F(Q) - F(Q ) + A.(Q-Q ) ; G(Q) - G(Q. .) + B.(Q - Q. .) 

i .j i ,J i iJ I *J 

*♦* + 

3F *> ^ 3G + 

A = — (Q. .) ; B = — (Q. .) 

3Q 1 ,J 3Q 




~ ■> 


(3.27) 
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where A and B are the Jacobian matrices evaluated at Q. . with the aid from 

i ,3 

Eq. (3.27), Eq. (3.26) can be written as 

(j 1 rv /v rv <v ->• 1 /v /v /v /v -f 4 1 

— Q. .+ (-y A + X 5 B)(Q. . - Q. . )+ (y A - x B) (Q . -Q . , > o 

dt i ,j ~ 5 i,j+l i,j-l „ ~ n n i+l,j i-l,j 

2^nJ 2 a$J 


which can be treated by Fourier analysis. Let 


(3.28) 


Q (t) = Q (t) = Qt (l 1) ( i 2) J (3.29) - ir<0 <ir, - ir<0 <tt (3.29) 

ij ij 1 2 

be the solution of Eq. (3.28) 

Substitution of Eq. (3.29) into Eq. (3.28) yields the eigenvalue problem 


i sin(e ) i*sin(e ) 

[si + L (y A - x B) + — (- y A + x B)] Q = 0 

A*J n n And 55 


(3.30) 


A nontrivial solution to Eq. (3.30) can be found only if 

y sin (e ) y sin(e ) - x sin(e ) x sin(e ) 

det [(-2 -1 i — )A + ( — L + _i L)b - (- isl) ] = 0 

L V fst / rv» /v y J 

A?J AnJ A£J AnJ 

(3.31) 

Since the original conservation law. Eq. (3.21) is hyperbolic, the matrices 

A and B by definition satisfy the condition 

det (oA + SB - xl) = 0 = > (x (ct,e)} n - all real (3.32) 

P P=1 

a, 8 real 

Thus, the eigenvalues s p ( E< 1* ( 3<31 ) are a11 purely imaginary and 
given by 

A 

s (0 , 0 ) = - i X ( a . 8 ) ; P = l,-««,n (3.33) 

p 1 2 p 


where 
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a - 


y sin (e ) y sin (e ) 
n l S L 


A?J 


AnJ 


x sin (e. ) x sin(e ) 
n i . 5 ^ 


A?J 


AnJ 


For the Euler equations the Jacobian matrices A and B are known analytically 
as well as the eigenvalues x (a,e). Thus, the local spectral radius at mesh 
point (i ,j) 

p i,j = maX l ^ S p ( 0 l 0 2 )l ^ 

P 9 9^ > ^2 


can be estimated by analytic means. The "local time step" scaling of Eq. 
(3.24) is, then defined by 


(// Jdgdn) p 
D 

ij 


— Q. 

i ,j dt l ,j 


+ 6 [- 
n 


(« y. .)(u F. . 
5 i.J n 


) + (5 x. .) (y 6 . .)] 
C i n i ;J 


+ <5 [<5 y Mil F. .) - (S x. .) (y 6 . J] = 0 (3.34) 

5 n i S J 5 1 S J n i S J C 1 »J 

which evidently scales the problem so that the local spectral radius is equal 
to 1 everywhere. It can be seen that this type of scaling does not affect a 
steady solution of the original scheme. 
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4. RESULTS AND DISCUSSIONS 

4.1 Sphere and Slender Body 

The inviscid flow over a sphere is chosen as the first application of the 
approach to the aerodynamic configurations. The simplicity of the sphere, 
both from grid generation point of view and flow field calculation, gives an 
advantage to the understanding of the approach. Since the approach has not 
been applied to a physical aerodynamics problem, the first application should 
be somewhat straight forward so that only difficulty lies in the treatment of 
the interface conditions. However, Euler equations model is not suitable for 
the high speed flow over sphere. This is because such a flow separates 
somewhere downstream and the Navier Stokes equations have to be used. In 

order for the Euler equations to be applicable to high speed flow, and yet 
simplicity of the configuration is maintained, a slender body is also 
considered. Solutions for flows over both configurations have been obtained 
at zero-degree angle of attack. Free stream Mach number for flow over sphere 
is 0.2 (nearly incompressible) while flow over the slender body has been 
investigated at a free stream Mach number of 1.5. The entire flowfield is 
divided into two subdomain at about the center of the configurations. Grids 
in both subdomains are of 0-0 type with different number of grid points, i.e., 
21x21x17 v.s. 33x33x17. Figure 4.1 and 4.2 illustrate grids used for sphere 
and slender body respectively. As mentioned earlier, the purpose for the 
investigations of flow over these configurations is to see whether the 
approach is feasible when applied to the simple CFD calculation. Solutions to 
these problems can be obtained by using single grids. The use of single grids 
is probably more efficient. Thus, solutions from single grid calculations are 
used as references to compare with solutions obtained by multiple grid 
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calculations. Wall pressure coefficient at the center line are compared in 
Fig. 4.3 and 4.4. Good agreement can be seen clearly from the comparisons. 
Even though, these comparisons are made for only simple cases but confidence 
in using the concept of multiple grids should be built up. Indeed, solutions 
to flow over these configurations can be obtained for various flow conditions 
and number of grid points. The same agreement is expected to be obtained. 
However, it is felt that problems with higher level of difficulty should be 
pursued. Thus, the use of different grid topologies for different subdomain 
grids is considered next. A Butler-wing described in the next section suits 
this purpose wel 1 . 


4.2 Butler Wing Configuration 

A Butler wing is a delta wing which was proposed by D.S. Butler [137]. 
The planform of the body is an isosceles triangle, and the leading edges of 
the wing lay along the Mach lines of the unperturbed stream. The first 20% of 
the wing is conical and the last 80% of the wing has elliptical cross section 
with increasing eccentricity along the x-axis. At the trailing edge, the 
elliptic cross section has infinite eccentricity and the last cross section is 
a straight line. Figure 4.5 shows a physical model of a Butler Wing. The 
semi major and minor axes are given by 


major axis (i.e., 

sime-span) = 

x/B 

0 < x < L 

(4.1a) 

minor axis (i .e. , 
center line) 

thickness on 
= x/e 


0 < x < 0.2 L 

(4.1b) 


■ - 6 [l - 

1 — » 

—1 

CVJ 

• 

O X 
CO 

1 • 

o 

X 

0.2L < x < L 

(4.1c) 

2 2 

where g = M - 1 
00 





Butler [137] 

has compared 

the experimental results for surface 

pressure 


with the theoretical results using the slender-body theory approximation to 
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simplify the inviscid equation of motion. Squire [138,139] has obtained 
experimental results for a Butler Wing with varying Mach number and angle of 
attack. Abolhassani et al . [140] have obtained numerical solutions on a 
Butler Wing by solving Navier-Stokes equations with the MacCormack time-split 
method. It should be mentioned that the experimental model is 120 mm. long 
and is constructed for M = 3.5. That is the semi-apex angle of the planform 
and of the initial conical nose is sin"*(l/3.5) = 16.602°. The model is 
mounted in the tunnel by attaching a string support of 12.5 mm. diameter to 
the lower surface. 

Even though solutions can be obtained using a single grid system, 
multiple grid solutions on a Butler Wing suit well for the purpose of this 
study. This is because grid of 0 type is suitable for the front part of the 
configuration while H-type grid seems to be a better choice near the trailing 
edge. The multiple grid system is composed of three subdomain grids as shown 
in Fig. 4.6. The H-type grid at the rear of the configuration is divided into 
two subdomain grids mainly to avoid programming difficulty. The Butler Wing 
configuration illustrates the usefulness of the multiple grids approach where 
grid topology is changed from one subdomain grid to another. 

Wall pressure coefficient at various flow conditions are plotted in Fig. 
4.3-4.10. In all plots, a solid line indicates wall pressure coefficient 
obtained from multiple grid calculations. Plots of wall pressure coefficient 
along the center line at 3.5 Mach number and zero degree angle of attack are 
shown in Fig. 4.3. In Fig. 4.7a, comparisons are made with ref. [136-139], 
Discrepencies near the nose region exist because grid is not fine enough to 
obtain the correct conical solutions there. Discrepencies at the rear of the 
wing are believed to occur because of the negligence of viscous effect in the 
Euler equations. Squire [138] has made similar conclusion about these 


discrepancies. Good agreement for the comparison between multiple grid 
solutions and solutions obtained from a single grid calculation is shown in 
Fig. 4.7b. This indicates that the use of multiple grids does not cause the 
discrepancies in Fig. 4.7a. 

The computed pressure coefficients for ten degrees angle of attack are 
plotted and compared with the Refs. 138 and 140 in Fig. 4.8. At 17%, 30%, 50% 
and 70%, chordwise positions, the pressure coefficients are plotted against 
the conical spanwise coordinate ^ tan Good agreement can be seen. On the 
thick sections near the nose, the pressure is highest on the centerline and 
falls toward the leading edge, whereas near the trailing edge the spanwise 
distribution is more 'wing like' with the maximum pressure at the leading 
edge. The changeover is shown by the pressure peaks in the pressure 
distributions at x/c * 0.5 and 0.7. Some discrepencies with the experimental 
results may be due to the fact that the lower surface of the experimental 
model is di stored to include a sting support. Results for flow at 2.5 Mach 
number are shown and compared with results from Ref. 138. Figure 4.9 shows 
results at zero degree angle of attack, whereas, results at ten degrees angle 
of attack is shown in Fig. 4.10. The same conclusion as in the case of 3.5 
Mach number can be made. 

The results obtained, thus far. demonstrate that the use of multiple grid 
approach is plausible and does not add significant error to the flow model 
equations even when grid topologies in subdomain grids are completely 
different. However, more complex problems should be investigated in order to 
be certain about the capabilities of the approach. This study may not yet 
indicate the usefulness or necessities of the multiple grid approach, since 
the construction of a single grid system can be made in all cases. In some 
applications, however, the construction of such a single grid system to cover 
the entire domain may not be possible at all. 
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5. CONCLUSIONS 

Solutions of Euler flow over aerodynamics configurations on multiple grid 
systems are presented. Some details on grid generation techniques and 
solution procedures are discussed. The concept of conservative treatment at 
the interfaces of various subdomain grid is also addressed. The solutions 
obtained from this study illustrate a promising future of multiple grids 
approach. It should be stressed, once more, that the main purpose of this 
study is to determine whether the use of multiple grids approach is feasible 
on these configurations, not to determine the characteristics of the flows. 
The use of multiple grid approach, however, should give the correct 
characteristics of such flows, as this study indicates. Thus, a single grid 
system can be constructed to solve for soultions on any of the configurations 
considered in this study. The solutions obtained from a single grid 
calculation can be used as references to compare with those obtained from a 
multiple grids calculations. This fact should not underestimate the 
usefulness of multiple grids approach. In some instances, the construction of 
a single grid system may not be possible at all. Even for simple 
configurations where a single grid system can be constructed, the use of 
multiple grids approach eliminates the difficulties that arise in the grid 
generation procedures. Also, experiences have shown that the use of multiple 
grid approach enhances the solution procedures. For example, the convergence 
to steady state of the solution to the equations of motion depends on many 
factors. One factor which is very important is the characteristics of a grid 
system, i.e. grid spacing, grid orthogonalities e.t.c. Good characteristics 
of a grid system is much easier to obtain in multiple grid approach as 
compared to a single grid approach. Experience from this study has indicated 
that efforts to be made to obtain a converged solution as a multiple grid 
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system is not as much as those to be made on a single grid system. Arguments 
can be made, however, that configurations used in this study are yet too 
simple to make any conclusion about this advantage of the multiple grid 
approach. More studies will have to be made to confirm it. 

The next step, which is under investigation, is to apply this approach to 
simulate the internal /external flow around a fighter-aircraft configuration 
(Fig. 5.1). This simulation should illustrate the usefulness of the multiple 
grid approach, since it becomes necessary to use different grid topologies 
between the internal and external flow. The findings from this study is 
expected to yield a significant contribution to the field of CFD. After the 
solutions of the Euler equations is sucessfully obtained, the applications to 
the Navier-Stokes equations, including chemical reacting and turbulent flows 
(providing that proper turbulent modelling is implemented), should be 
possible. 
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(a) Physical domain 


(b) Computational domain 


Fig. 2.2 Physical versus Computational domain. 












Fig. 2.7 An 0-0 mapping for a wing-fuselage configuration. 













Fig. 3.2 A solid wall boundary. 
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Fig. 3.4 Conditions at the inflow /outflow boundaries. 
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Fig. 3.6 A typical patched interface. 
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(a) Grid 1 



(b) Grid 2 



(c) Symmetry plane 


Fig. 4.1 Grid system for a sphere 
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(a) Grid 1 



(b) Grid 2 

Interface 
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(c) Symmetry plane 


Fig. 4.2 Grid system for a slender body. 
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Fig. 4.3 The pressure coefficient on the surface of a sphere. 
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Fig. 4.4 The pressure coefficient on the surface of a slender body. 
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(a) O-Grid 


(b) H-Grid 



(c) Symmetry plane 


Fig. 4.6 Grid system used for a Butler-Wing, 


Present Butler wing 
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Fig . 4.7(b) Comparison of C p with the single grid calculation. 
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Fig . 4.8(a) C p on the surface of a Butler-Wing at different cross cut (M = 3.5, 
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Squire (1981) M„= 3.5 
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Fig. 4.8 (d) x/c = 70% 
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Chord Length 

Fig. 4.9 C p on the surface of a Butler- wing(M = 2.5, a = 0°) 
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Fig. 4.10(a) C p on the surface of a Butler-Wing at different cross cut 
(M = 2.5, a = 10°) 
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Fig. 4.10(b) x/c = 30% 
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Fig. 4.10(d) x/c = 70 % 



